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ABSTRACT 

We report on the development of Mezcal-SRHD, a new adaptive mesh refinement, special relativistic hydro- 
dynamics (SRHD) code, developed with the aim of studying the highly relativistic flows in Gamma-Ray Burst 
sources. The SRHD equations are solved using finite volume conservative solvers, with second order interpola- 
tion in space and time. The correct implementation of the algorithms is verified by one-dimensional (ID) shock 
tube and multidimensional tests. The code is then applied to study the propagation of ID spherical impulsive 
blast waves expanding in a stratified medium with p oc r"* , bridging between the relativistic and Newtonian 
phases (which are described by the Blandford-McKee and Sedov-Taylor self-similar solutions, respectively), 
as well as to a two-dimensional (2D) cylindrically symmetric impulsive jet propagating in a constant density 
medium. It is shown that the deceleration to non-relativistic speeds in one-dimension occurs on scales signif- 
icantly larger than the Sedov length. This transition is further delayed with respect to the Sedov length as the 
degree of stratification of the ambient medium is increased. This result, together with the scaling of position, 
Lorentz factor and the shock velocity as a function of time and shock radius, is explained here using a simple 
analytical model based on energy conservation. The method used for calculating the afterglow radiation by 
post-processing the results of the simulations is described in detail. The light curves computed using the results 
of ID numerical simulations during the relativistic stage correctly reproduce those calculated assuming the 
self-similar Blandford-McKee solution for the evolution of the flow. The jet dynamics from our 2D simulations 
and the resulting afterglow lightcurves, including the jet break, are in good agreement with those presented in 
previous works. Finally, we show how the details of the dynamics critically depend on properly resolving the 
structure of the relativistic flow. 

Subject headings: gamma rays: bursts - hydrodynamics - methods; numerical - relativity 



1. INTRODUCTION 

Gamma-Ray Bursts (GRBs) are the most electromagnet- 
ically luminous explosions in the Universe. Their non- 
thermal and highly variable gamma-ray emission implies 
that the emitting region must be ultra-relativistic - mov- 
ing with a very large Lorentz factor, typically > 100 and 
sometimes as high as > 10^^, in or der to avoid excessive 
pai r production at the source_ (e.g.. iLithwick & Sari l 120011: 
Granot. Cohen-Tanugi_&do Couto e Silv a 2008; Abd o et al] 
2009allbHAckermann et all 120101) . At suflicientlv large dis- 
tances from the source the GRB outflow decelerates as it 
drives a strong relativist ic shock into t h e surrounding medium 
(for reviews see, e.g., IPiraiil l2005al: iGranotI l2007l) . Syn- 
chrotron emission from this long lived external shock powers 
the GRB afterglow, which is observed in the X-rays, optical or 
radio, typically over days to months after the prompt gamma- 
ray emission. The peak frequency of the afterglow emission 
shifts to lower energies as the afterglow shock decelerates by 
sweeping up the external medium (Rees & Meszaros 1992|). 
This picture of a decelerating relativistic expansion of the 
emitting region during the afterglow phase is supported by 
direct measurements of the afterglow image size at late times 
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in the radio, using very long base-line in terferometric tech- 
niques, for GRB 030329 a t z = 0.1685 dTavloret al.l 120041 
l200llPihis^om et al.ll2007l) . 

GRB activity manif ests itself over a dynamical range o f 
~ 13 decades in radius (iGehrels. Ra mirez-Ruiz & Foxll2009l) . 
The phenomena involves different stages, which are usually 
modeled separately because of their complexity. Let us con- 
sider these stages in turn, working from the small scales to the 
large scales. 

1.1. Jet Production and the Central Engine 

GRBs divide into t wo classes according to t heir duration 
and spectral hardness (iKouveliotou et alJll993h . Long dura- 
tion GRBs (lasting > 2 s) are associated with Typ e Ic core col- 
lapse SNe, and thus to the d eath of massive stars dStanek et alj 
I2OOI iHiorth et al.l I200I IWooslev & BioomI l2006h . while 
the nature of short dur ation GRB (lasting ^ ^ Pro- 
genito rs is still debated jLee & Ramirez-RuizI l2007t iNakaij 
l2007h . the most popula r model involving the binary merger 
of two compact stars (iPaczynski ll986tlEichler et aP 1 19891 : 
iNaravan. Paczvnski. & Piranlll992 ). 

In the collapsar model for long GRBs ( IWooslevll993h . dur- 
ing the collapse of a massive Wolf-Rayet progenitor star a 
black hole is formed, which rapidly accretes stellar envelope 
material, launching a relativisti c jet that penetrates the star and 
eventually powers the GRB dRamirez-Ruiz. Celotti & ReesI 
l2002h . It has been modeled using numerical simulations, 
where a jet is usually injected as an inner boundary condition 
at the center of a collapsing massive star, and bores its way 
out of the progenitor star's envelope (MacFa dyen & Woosleyl 
Il999t IZhang. Woosley & MacFadyen 2003) . Some simula- 
tions include a magnetic field (in an ideal magnetohydrody- 
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namical fra mework) and recently ad ded a general relativistic 
framework (iMizuno et al.l l2004allbl: lHa\ylev & Krolig 120061: 



McKin nev' '2006^, 'Nagataki et al.n2007l; iTchekhovskov et al.l 
20081 [Nagataki 2009;^ Barkov & Baushevll201 ll) . The alter- 



native model for the central engine of long GRBs featur- 
ing the formation of a millisecond magnetar (i.e. a very 
rapidly rotating highly magnetized ne utron star: lUsovll 19921) 
has a l so been studied numerically (" Komissarov & Barkovl 
120071: iBucciantini et al.l 12007. 2008, 2003). Binary merger 
simulations of two neutron stars or a neutron star and a 
black hole were performed in the context of short GRBs 
(iLee & Ram irez-Ruiz 2002: Rosswog & Ramirez-Ruiz 2002; 
Rosswod 12005; .Faber et al., ,2006^: ,Oechslin & Janka ,20061: 



Rezzolla et al. Recent general relativistic magneto- 



hydrodynamics (MHD) simulations show that a relativistic 
jet can naturally fo rm in such a scenari o, which may indeed 
power short GRBs (iRezzoUa et al.l20ril) . Similar simulations 
of relativistic jet formation from accretion onto a black hole 
are routinely performed also i n the contex t of active galactic 
nucle i or micro-quasars (e.g.. lMeieill2003l:lKrolik & HawlevI 
120101) . Different processes have been suggested to accelerate 
and collimate the jet: (i) thermal energy injected into the jet 
by annihilation of neutrinos and anti-neutrinos from an accre- 
tion disk (e.g.. iFrver & WooslevHl998t iPopham et al.lll999 



Rosswog, Ramirez-Ruiz & Davies 2003: Lee et al. 12004 : 
Lee & Ramirez-Ruiz 2006: Chen & Beloborodov 20()7|); 00 
rotational energy extracted from the central black hole 
through the Blandford-Znajek effect (iBlandford & Znaiekl 



glow stage are usually done separately from the earlier stages 
(of the jet formation, acceleration and collimation), in order to 
simplify these challenging numerical computations, which in- 
volve a very large dynamical range. The most common initial 
conditions for simulations of the GRB jet during the afterglow 
stage are a conical wedge of half-opening angle 6o taken out of 
the spherical BMK solution (though in some cases a relativis- 
tic cold shell or blob is used instead). Since the angular size 
of regions that are casually connected in the lateral direction 
is ~ 1 /r, such a BMK wedge should not evolve significantly 
while its Lorentz factor is F » \ suggesting that the sub- 
sequent evolution should be insensitive to the exact choice of 
initial Lorentz factor Fq in the limit where Fq » 0q ' . 

For an ultra-relativistic blast wave most of the energy in the 
shocked (downstream) region is within a thin layer behind the 
shock transition, whose width is A ~ Q.IR/F^ in the lab frame 
(i.e. the rest frame of the external or upstream medium, which 
in our case is also that of the central source), which is hard to 
resolve prope rly for large initial Lorentz factors (see, e.g., 
lGranotll2007l) . Therefore, most simulations use Fq^o ~ 3 - 4 
rather than the ideal choice of Fq^q ^ 1, along with values 
of 6() that are not very small (usually 6*0 = 0.2 and Fq ~ 20), 
despite the actual initial Lorentz factors at the onset of the 
afterglow are estimate d to be at least a few hundred (e.g., 
iLithwick & Sarill2001l) . while the values of inferred from 
afterglow observations (e.g., [Frail et al. 200 li) can be as low 
as ~ 0.03 - 0.05 (or a high as > 0.5). 



through the Blandiord-ZnaieK etiect OBlandiord & ZnaieKI 
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(iii) rotational energy extracted from the accretion disk, 
coupled with a dynamica lly important magnetic field 
jBlandford & Payne 1982: ProgaetalJ 120031: iLvnden-Belll 
[2003CiUzdenskv&MacFadveft200"6l) . 

1.2. Jet Expansion and Deceleration 

Once the GRB outflow transfers most of its energy to 
the shocked external medium it becomes dynamically sub- 
dominant and the flow becomes insensitive to the ex- 
act composition or initial radial structure of the original 
outflow. A t this stage a spherical flow approaches the 
lBlandford & McKee (1976) self-similar solution (hereafter 
BMK), losing memory of the initial conditions and retain- 
ing memory only of the total energy. The complete evolution 
of a spherical relativistic fireball, including the acceleration, 
coasting and deceleration phases, has been studied numeri- 
cally by Kobavashi, Piran & Sari ( 1999) by using one dimen- 
sional (ID) spherical simulations. 

When a non-spherical relativistic outflow (or jet) decel- 
erates, to zeroth order it locally resembles a section of the 
spherical BMK solution characterized by the local value of 
the energy per solid angle or isotropic equivalent kinetic en- 
ergy, iiiciso- Once the Lorentz factor F drops to 0y', where 
Oq is the initial half-opening angle of an initially uniform 
jet with sharp edges, the jet becomes causally connected 
in the lateral direction and can in principal start spreading 
sideways significantly. Simple analytic mo dels argue that 
it should indeed quickly spread sideways (iRhoadsl ll997L 
119991: ISari. Piran & Halpernll 19991) . while numerical simula- 
tions show that the lateral spreading is much more modest, 
and the flow retains memory of 0q for a long time, which for 
typical values of Oq in GRBs lasts u p to the non-relat ivistic 
transition time (Granot et al."2001': Granot & K 6niglll2003l: 
ICannizzo et"a!1[2004; .Zhang & MacFadven,2009,) . 

The numerical simulations of jet dynamics during the after- 



Since afterglow emission is thought to be predominantly 
synchrotron radiation from the shocked external medium, 
then accurately inferring the properties of the original rel- 
ativistic outflow and the external medium from afterglow 
observations requires an accurate modeling of the dynam- 
ics. The jet numerical simulations and calculations of the 
corresponding afterglow emission (Granot et al . 2001) have 
recently been extended to well wi t hin the non-relativistic 
stage (e.g.. IZhang & MacFadveij 120091: Ivan Eerten et alj 
I2010t IWvgoda et aljboill: Ivan Eerten & MacFadvenll201 It) . 
Following the dynamics from a highly ultra-relativistic initial 
Lorentz factor (Fq > 20, for which Ao//?o ~ 10-4(Fo/30)-2) 
down to highly Newtonian velocities (v < 0.01c) requires a 
very large range of spatial scales, for which an adaptive mesh 
refinement (AMR) code is necessary in order to properly cal- 
culate the multi-dimensional flow dynamics. iGranot et alj 
(2001) were the first to study this problem numerically by 
using multi-dimension numerical simulations and found that 
the GRB jet sideways expansion is slower than expected 
fro m analytical models. These results were later confirmed 
by IZhang & MacFadvenI (|2009|) . who followed the evolu- 
tion of the GRB jet up to the non-relativistic phase by run- 
ning high resolution two-dimensional (2D) simulations. Sim- 
ulations using similar ini tial conditions were also run by 
iMeliani & KeppensI (1201 Ol) who found that the shock front be- 
comes unstable at high values of the Lorentz factor, F > 15, 
but the instabilities quickly decay when the jet decelerate to 
F< 10. 

All the multi-dimensional numerical simulations of af- 
terglow jets have so far assumed a uniform external 
medium, even though a stratified external medium is ex- 
pected fo r the stellar win d of a massive star long GRB pro- 
genitor ([C hevalier & Li 2000| iPanaitescu & Kumad [2OO0I: 
iRamirez^Ruiz et al.. 2001, .20051) . This was partly motivated 
by the faster deceleration of the afterglow shock with ra- 
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dius in a uniform external medium compared to a strati- 
fied one, which reduces the required dynamical range of 
the simulations. Moreover, magnetic fields may also aff'ect 
the jet dynamics (in additi on to their effect on the after- 
glow synchrotron radiation). iMimica et al.l (120091 ^20 10) have 
used ID simulations to study the deceleration of magnetized 
GRB ejecta propagating into a uniform ambient medium, 
and showed that while the late evolution of strongly magne- 
tized shells resembles that of hydrodynamic shells, the mag- 
netization plays an im portant role into the onset of the for- 
ward shock emission. iMimica & GianniosI (1201 lb computed 
the afterglow emission produced by a GRB ejecta deceler- 
ating into a realistic external medium by running ID spher- 
ical simulations. However, multi-dimensional simulations 
are necessary in order to fully capture the magnetic field 
dynamics, as for instance the generation of turbulence by 
the m agnetohydrodynamics Kelvin -Helmholtz (|Z hang et al. 
2009) or Richtmyer-Meshkov dGoodman & MacFadvenI 
20081) instabilities, and the co nsequent mag netic field am- 
plific ation ("Inoue et al.' 20 11; Mizuno et al. '2011'). Actu- 
ally, iGranot. Komissarov & Spitkovsky (2011) have recently 
shown that even in ID one cannot realistically model the de- 
celeration stage separately from the acceleration stage if the 
outflow is initially highly magnetized and accelerates under 
its own magnetic pressure. Instead, a full simulation of the 
acceleration and deceleration is needed, requiring a very large 
dynamical range that is numerically challenging. 

With the aim of addressing these questions, and perhaps 
also possible applicability to earlier stages of the jet dynamics 
(such as its acceleration or propagation within the progenitor 
star), we have developed a new AMR, relativistic hydrody- 
namic code. While the code developed is sim ilar in several 
aspec ts to previous SRHD-A MR code s (e.g., iHughes et alJ 

^advenl 



2002t lAnninos et a .1 120051: Zhang & MacFadvenI 120091 
Meliani et al.1 120071: iMorsonv et al.ii2007l: IWang et al.ll200 
we consider it important to present a detailed, self-contained 
description of the hydrodynamic code as well as the match- 
ing radiation code, along with detailed tests. The paper is 
organized as follows. ^ and ^ describe in detail, respec- 
tively, the SRHD code and the radiation code used to calculate 
the observed afterglow emission (by post-processing the out- 
come of the SRHD simulation). Standard tests used to verify 
the SRHD code are presented in the Appendix, while the cor- 
rect implementation of the radiation code is discussed in ^ 
^presents a detailed study of the propagation of a relativis- 
tic, purely hydrodynamic ejecta into a one-dimensional strat- 
ified medium as well in a multi-dimensional homogeneous 
medium together with the resulting lightcurves. Finally, ^ 
presents our conclusions. Simulations of the propagation of 
jets into a stratified medium and the inclusion of magnetized 
flows will be addressed in future work. 

2. NUMERICAL CODE 

2.1. Relativistic Hydrodynamics equations 

The special relativistic hydrodynamics (SRHD) equations 
in conservative form (e.g.. iAnile..l989,) can be written as fol- 
lows: 



f +V.(Dvt).0 

|.v.(^..,x).o 

+ V ■{tv + pv) = 

ot 



(1) 
(2) 
(3) 



where p is the thermal pressure, v = /?c is the flow velocity (c 
being the speed of light), and I is the identity matrix. These 
equations represent the conservation of rest mass ([TJ, momen- 
tum (l2]i, and energy The conserved variables {D, 5*, t) 
correspond to the lab frame rest mass, momentum, and energy 
(excluding rest mass) densities, respectively. They are related 
to the primitive variables (p, v, p) by the following relations: 



D = pr, 

§ = DhTv , 



DhYc^ -p- Dc^ 



(4) 
(5) 
(6) 

where F = (1 - ffy^^^ is the Lorentz factor, p is the proper 
rest mass density, and h is the specific enthalpy. The SRHD 
system of equations is closed by the equation of state, relat- 
ing h to p and p. Note that by explicitly subtracting the rest 
mass in the definition of the lab frame energy density t in 
equation (|6]l, the non-relativistic hydrodynamic equations are 
properly recovered when /? <*: 1 . 

2.2. Integration methods 

The SRHD equations ([Tli-(l3]l form an hyperbolic system of 
equations and can be solved by using methods similar to those 
developed for cl assical non -relativistic gas dynamics (for a 
review see, e.g., lTordl200"8h . Without loss of generality, the 
solution of the hyperbolic system of equations 

du -> 

^-.V./ = 0, (7) 

is given in ID (the generalization to multi-dimensions is 
straightforward) by: 



Ax/ 



-,;i+l/2 



pn+l/2x 
■ ^/-l/2 ■ 



(8) 



where x, represent the position of the center of the cell ; with 
volume Ax; — x,+i/2 - x,_i/2, x,±i/2 are the positions of the 
interfaces between the cells x, and x,±i, and 



, x)dx , 



1 nxi^ii2 
« = — «,(f", 

Ax,- 

1/2^ = ^J^ fit,Xi^ll2)dt , 



(9) 



(10) 



are the volume average of the conservative variables and their 
time-averaged fluxes. 

While equation ([8]l represents an exact solution of the 
corresponding partial differential equation, an approxima- 
tion is introduced when the fluxes (equation [TOb are com- 
puted. Because an exact solver is in general very expen- 
sive, in the current version of the code we have implemented 
the simple and comput ationally efficient rel ativistic exten sion 
(ISchneider et alJfT993l) of the HLL method (lHartenlll983l) . 

It is well-known that the HLL method does not resolve 
properly the contact discontinuity, and it has an intrinsic high 
level of numerical diffusivity, while for instance other meth- 
ods (e.g. the HLLC method, Mignone & Bodo 2005) prop- 
erly reconstructs the contact discontinuity, producing results 
with significantly lower dissipation. On the other hand, being 
more diffusive, the HLL method is also more "robust", very 
rarely producing unphysical results such as negative pressures 
or imaginary Lorentz factors. In addition, a low dissipation 
method may produce undesirable effects, such as a "carbun- 
cle" artifact along th e axis of propaga tion of strong shocks 
(see the discussion bv lWang et al.ll20()8l) . 
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Second order accuracy in time and space are obtained by 
employing a Runge-Kutta integrat or and by a spa tial recon- 
struction of the primitive variables dvan Leeij|1979h . except in 
shocks where the methods drops to first order (in space) by 
a limiter. Different limiters are implemented, including the 
"minmod" (being the most diffusive), UMIST, Superbee and 
the less diffusive "monotonized central difference" limiter 

2.3. Extension to cylindrical and spherical coordinates 

The extension to cylindrical and spherical coordinates is 
treated very carefully in the code. For instance, in two- 
dimensional (r, ff) spherical coordinates, the equations read: 



dU _^ 1 d{r-F) 



1 d{Gsme) S 



(11) 



dt r^ dr rsinO dO r 

where U, F, G, S can be easily derived from equations ([T]i- 
([3]). We note that a simple cell-center discretization of this 
system of equations introduces large numerical errors when 
differencing. In particular, it does not preserve stationary ini- 
tial conditions to machine accuracy. As an example, if one 
assumes static initial conditions, as dp/dr = 0, p constant 
and V = 0, these are preserved in the code if, e.g., the rela- 
tion (easily derived from the 0-component of the momentum 
equation) 

1 dipsinO) pcosO 

' - ' (12) 



r sin 6 dO 



r sinO ' 



is held to machine accuracy. A simple centered discretization 
gives 

1 sin 0j+ 1 /2 - sin 1 /2 cos Oj 



sm( 



sin 6: 



(13) 



where Oj is evaluated at the center of the cell, while 0y±i/2 at 
the interface between different cells, and it does not preserve 
the initial conditions. 

A way to minimize numerical errors when differencing 
equation ( fTTT ). especially near coo rdinate singularities, is by a 
finite volume discretization (e.g.. lFalleiri99U iLi & Lill2003h . 
that is by averaging the variables over the cell volume. Given 
for instance the cell centered in (/, j) and with nodes located 
at (/ ± 1/2, j + 1 /2), the value of of the quantity A averaged 
over the cell volume is given by 



(A) 



Jsinddd jAr-dr 
JsinOde Jr^dr 



(14) 



With this definition, radial and polar derivatives are approx- 
imated by (taking A = 77^^ and A = 7^17^^^^^ respec- 
tively): 



1 d{r^F) diir'-F) 



1 



r2 dr 6i(ry3) 
d{Gsin9) 6{G sin 6) 6{r^/ 2) 



rsinO 



86 



(15) 



-S{cos0) 6(r^l2>) 

where (5,(/) = fi+iji - /i-1/2, while the source terms are dis 
cretized by assuming (taking A - ^ and A - 
tively): 



sin^ 



re spec - 



1 



Si{r-/2) 



cos 6 6{sin ff) 



(16) 



r 5,(rV3) ' sine -6(cos0) ' 

It is easy to verify that, written in this form, equation fT% 
preserves static initial conditions to machine accuracy. 



2.4. Equation of state 

The equation of state relates the enthalpy to the pressure 
and d ensity. In the case of a relativistic perfect gas it takes the 
form (lSvngdll971l) 



^2(1/0) 



(17) 



where = p/(pc^), and Ki is the /*-order of the modified 
Bessel functions of the second kind. 

As the evaluation of the enthalp y from equation ([17]) is com - 
putationally expensive (see e.g. iFalle & Komissarovlll996h . 
simplified relations have been used, the simplest being the y- 
law equation of state (EOS) 



h = 1 



1 



-0. 



(18) 



with a constant value of the adiabatic index y fixed and equal 
to 4/3 or 5/3, valid only in the limit of ultra-relativistic or 
s ub-relativistic fluid s, respectively. 

iMignone & Bodol (i2005.) proposed the EOS (see also 
lMathewslll97°li) 



(19) 



which in addition to approximating equation ([TtI i within 2%, 
also satisfies thejaub (1948i) inequality 



{h - 0)(/i - 40) < 1 , 



(20) 



in accordance wit h relativist i c kinet ic theory. 

More recently, iRyu et all (|20 06b proposed a simpler and 
better approximation to the Synge EOS (accurate to within 
0.5%), which also satisfies the Taub inequality (equationl20ll. 
given by 

^^^60^^ (21) 
30 -F 2 

The implementation of these EOS is straightforward, and un- 
less stated otherwise, in this paper we use the one derived by 
iRyu et alJ (l2006h . 

2.5. Converting conserved to primitive variables 

The increased level of complexity in solving the SRHD 
equations when compared to the corresponding non- 
relativistic hydrodynamics equations arises mainly from the 
lack of simple closed expressions relating conserved (t, 5*, D) 
and primitive {p, v, p) variables. This requires the primitive 
variables to be computed from the conserved variables by a 
non-linear iteration; 

Among others, iNoble et"aI1 (l2006h studied several algo- 
rithms to convert conserved t o prim itive variables for the case 
of a y-law EOS. iRvu et al.l ( I2OO6I) . for the EOS defined in 
equation (I2TI 1. applied a Newton-R aphson method to an 8- 
th ord er equation dependent on F. iMignone & McKinneyl 
( 120071) . for the case of relativistic MHD with a general equa- 
tion of state, derived an equation for W - Dph, and evaluate 
W by a Newton-Raphson iterative scheme, with the deriva- 
tive dW/dp given by using thermodynamics relations. Here, 
we present a different implementation. Taking advantage of 
the existence of a relation between the specific enthalpy /; and 
- plipc^), we solve the system of equations (|4]|6l) as a func- 
tion of by using a standard Newton-Raphson method, and 
we then determine the other variables. 
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First, squaring the momentum equation {S k - DhYvk) we 
get: 

(22) 

with h - h(@). From the definition of specific enthalpy it 
follows that h > I. Therefore, equation ( l22l i leads to the fol- 
lowing inequality (e.g. Schneider et al. 1993) 



(23) 



By using the relation p - DQc^/F, we can then derive from 
the definition of energy density (excluding rest mass, i.e. t - 
DhTc^ - p - Dc~) the following identity 



/(0) = /i(©)r(©)- 



© 

n©) 



(24) 



Equations ( |22] | and ( l24l i are then used, together with a stan- 
dard Newton-Raphson method, to determine ©, with df/d@ 
given by 



df(@) ^h'l © r^-i \ 1 
d@ r\ /j r2 / r 



(25) 



where the relation F' = -h'{r^ - l)/{hr) has been used (de- 
rived from equation l22Ti. and h' = dh/d&. The derivative 
dh/d& depends on the particular EOS used, and can be de- 
termined both an alytically or numerically. In the case of the 
lRvuetal.l(l2006l) EOS (equation|2B, /i' = 4 - 6/(3© + 2f. 

We also note that df(&)/d@ > for every value of © (for 
the EOS considered here). Therefore, as /(© ^ oo) > 0, a 
solution for the equation /(©) = exists if /(© = 0) < 0, 
which implies the relation 



(26) 



must hold in order to allow a solution with physically accept- 
able values of F and p (that is, real values of F > 1 and p > 0). 

As we have shown, this method can be easily applied to any 
equation of state of the form h - h{&). Furthermore, the guess 
used by the Newton-Raphson method (NRM) is provided by 
simply assuming © = 0. In this case, setting a tolerance of 
10""^ into the Newton-Raphson solver, the method converges 
typically within ~ 5 iterations. In very rare cases when the 
NRM fails to converge, a bisection method is used instead. 

2.6. Adaptive mesh refinement 

We have implemented the SRHD equations in the frame- 
work of the adaptive mesh refinement code Mezcal. In the 
code, a basic Cartesian grid is built at the beginning of the 
simulation, and it is refined based on the initial conditions and 
the subsequent evolution of the flow. The uniform version of 
the code has been used in the past to simulate MHD jets (e.g., 
iDe Colle & Ragalllool l2006tlDe Colle et al.llIOOSh . 

In the Mezcal code, the computational grid is divided in 
"octs" (or blocks) of 2"''™ cells, where «dim is the number of 
dimensions of the problem. Each block has a series of point- 
ers to its vertexes, and each vertex has pointers to the octs 
sharing that particular vertex. In this way, neighbor octs (both 
along the axes and the diagonal direction) can be easily lo- 
cated in the grid, facilitating the computation of the MHD 
solver (that, in staggered mesh methods, is based on determin- 
ing electric fields at the cell vertexes). At a given time, each 
position on the grid is covered by only one cell, i.e. there are 
no pointers between "parent" and "sibling" usually present in 



other tree-AMR codes (e.g. iBerger & OhgedlI984t iKhokhlovl 
Il998h . Furthermore, there are no ghost cells in any of the 
blocks. Although the use of pointers causes a small mem- 
ory overload (corresponding to 4 integers per cell in 3 dimen- 
sions), that is largely compensated by the fact that, due to the 
small block size, the grid covers only regions that effectively 
need to be refined. 

At every timestep, all blocks are swept, and they are re- 
fined/coarsened if a user defined criterion is fulfilled. Typi- 
cally, this criterion is based on the first or second derivative of 
some variable, but more complex criteria can be easily imple- 
mented. Once a list of blocks flagged for refinement has been 
formed, the grid is checked for consistency. As the code main- 
tains a maximum ratio of 2 in the size of neighbor cells, all 
coarser neighbors of blocks are flagged for refinement. When 
a block is refined, 2"*"' new blocks are created, and the parent 
block is eliminated. To avoid excessive memory fragmenta- 
tion, the block lists are periodically reordered. 

Coarsening is allowed only when the 2"'"" neighbor blocks 
(previously produced by refining the same parent block) are 
marked for derefinement during the same timestep. We use 
zeroth-order interpolation when refining, and we integrate the 
conserved variables over the v olume when co arsening, fol- 
lowing the strategy presented bv lLi & O(l2003l) . 

To evolve the hyperbolic equations, the code employs a 
timestep common to all grid levels. While the use of a global 
timestep may potentially produce an important computational 
overload (as large of 5 0%, depending on the problem, see e.g. 
iDursi & Zingalell2003h with respect to using a local timestep, 
the local time step method can represent an important bottle- 
neck for parallelization, as blocks on different levels must to 
be evolved sequentially (and not in parallel). The fluxes are 
computed by locating the neighbor blocks, and considering 
the cells sharing the same faces. When two blocks with dif- 
ferent levels of refinement share the same face, 2"'"°'"' fluxes 
are computed between the 2"'""'"' cells located on the higher 
level block and the cell part of the block at the lower level of 
refinement. The fluxes are then added to the conserved vari- 
ables of the cells sharing the common boundary. 

The Mezcal code is parallelized by using the Message Pass- 
ing Interphase (MPI) library. The communication time is 
minimized by scheduling it in parallel with the calculation of 
the fluxes. This is done by first computing the fluxes between 
blocks located in each process, and then, once the commu- 
nication phase is completed, computing the rest of the fluxes 
(between blocks "inside" each process and ghost blocks). The 
load balanci ng is achieve d by ordering the blocks by a space- 
filling curve (ISaganI 19941) . dividing the total number of blocks 
between the diff'erent processes, and moving blocks between 
unbalanced processe s. In the code , the Morton and the Hilbert 
space-filling curves (ISaganlll994h are implemented. The load 
balancing is typically applied every ~10 timestep, and repre- 
sents an overload of ~ 1 % of the total computational time. The 
parallel scaling of the AMR code is under evaluation and will 
be presented elsewhere. 

3. CALCULATION OF THE EMITTED RADIATION FROM A 
HYDRODYNAMIC SIMULATION 

3.1. Calculation of the observed fiux density 

Here we provide a detailed derivation of the procedure 
required to calculate the radiation emitted from a r elativis- 
tic source, following iGranot & Ramirez-Ruiz (2010), which 
is based on previous work (iGranot. Piran & Sarii il999aiibl : 
IGranot & Koniglll2003l: iKumar & Granotll2003h . 
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where t~ is the coordinate time at the source's cosmological 
frame. 



Fig. 1 . — The contribution of a volume element dV to the flux observed by a 
distant observer is dFy(hi) = Iv(n) cos 6,^ rffisd ~ Iv(n)dQ.si, where 6si is the 
angle between the direction opposite to that at which the detector is pointing 
{hj = Z in the figure) and the local direction from a small emitting region 
within the source (of volume dV) to the detector. Since the observer is far 
away, the direction of emission in the observer frame is almost parallel to the 



The geometry of the problem is shown in Figure [T] We de- 
note with 0sd the angle subtended by the direction «d of the ob- 
server (perpendicular to the differential area dA at the detec- 
tor, and opposite to the direction at which the detector is point- 
ing) and the local direction h from the relevant (contributing) 
part of the source to the observer. In practice almost always 
flsd 1, as the source size is much smaller than the distance 
from the source to the observer, so that cos 0sd ~ 1 ■ We also 
define (iOsd = d(p^id cos 0sd as the differential solid angle sub- 
tended by the contributing portion of the source as viewed by 
the observer Our aim is to calculate the observed flux density, 
Fy - dE /dAdvdt, which is the energy per unit area, frequency 
and time in the direction hd normal to dA. From the definition 
of the angular distance to the source, dAiz), where z is the 
cosmological redshift, we have c/Qsd = dS ±/d^, where dS ± 
is the differential area in the plane of the sky (normal to «) 
sustained by the source. The angular distance to the source is 
related to the luminosity distance: diiz), hy d^ - (1 + zY^di. 

The differential contribution to the flux can be written as 
dFy{hd) - /v(«) cos 0sd <^i^sd ~ Iv(n)dQ.^i - lydS ±/d\. Here 
Iy(n) = dE I dAdQ.dvdt is the specific intensity (the energy per 
unit area, time and frequency of radiation directed within a 
small solid angle c/O, which is centered on the direction «), 
and should be evaluated at the location of the observer. 

For an optically thin source = J jvds,, where yV. = 
dE~l dV,dQ.,dv^dt~ is the emitted energy per unit volume, 
solid angle, frequency and time, while ds. is the differen- 
tial path length along the trajectory of a photon that reaches 
the observer at the time fobs when Fy is measured (the sub- 
script z here denotes quantities measured in the cosmolog- 
ical frame of the so urce). Since ly/v^, jy/v^ and ds/v 
are Lorentz invariant (iRvbicki & Lightmanlll979l) . we have 
ly - (v/v.)^Iy^ = (1 + zy^ J jy.ds-. Therefore, dFyih^) - 

IydS±/d^ = jyMV^il +z)ld\, where dV~ - dS ±dsz is the 
volume element in the source cosmological frame. Here 

jy_ = [r(l - n ■ li)Y^ fy, is measured in the source (cosmo- 
logical) frame, while f, is measured in the (comoving) rest 

frame of the emitting material, which expands at a velocity 0c 
in the source frame. Altogether, this gives^ 

T7 U -^ (1 +^) r ^4 A, n ■ r Jobs \ 

Fyitobs,n) ^ dxSlt, 

dHz) J \ c l+z 



v' = {]_+z)r(l-h-/3)v, 



fobs = (l+z)|fc- ^ , (28) 



and fobs = corresponds to a photon emitted at the origin 
(r = 0) at fj = 0. Since d'^x - dt-dV- - dt~dS±ds. - 
dUdS^ds'(vJV) = dudV IT(\-h-p) &nAAnj[,dV' = dL'^, = 
An{dE' I dil' dv' dt') is the differential of the isotropic equiva- 
lent spectral luminosity in the comoving frame, equation dZTI i 
can be rewritten as 

Fy{tobs, n) = -—^ \dh5\ U X 

I 1 p l+z 



Andliz) ■ 



■I 



dL' 



r3(l -n-yS)3 



(29) 



There are two main approaches to calculate Fy from the re- 
sults of a numerical simulation. The first one relies on numer- 
ically calculating along different lines of sight (i.e. trajec- 
tories or world lines of photons that reach the observer), and 
then comp uting dFy - lydS ±/dz. This was applied both in 
analytical (iGranot. Piran & Sarilll99 9b: Granot & Sari 200^ 
and in nurn erical (ISalmonson et al.i i2006t: ivan Eerten et al] 
|2010l 120111) calculations. Its main advantages are that it can 
properly handle the optically thick regime, where the radia- 
tive transfer equation is solved (analytically or numerically) 
along each line of sight, and that it provides the observed 
image of the source (i.e. ly on the plane of the sky) as a 
by-product, since it is used when calculating Fy. Its main 
disadvantage for numerical simulations is that it requires ac- 
cessing many different "snapshots" of the simulation results, 
corresponding to different lab frame times t^., for calculating 
each value of ly, as it requires integration along the trajecto- 
ries (or world-lines) of photons that reach the observer. The 
second approach, we adopt here, avoids this d ifficulty, and 
was already used in several previous studies (iGranot et alj 
20011 12OO2L IGranot & Koniglll2003t iKumar & Granoll 120031: 



Nakar & GranotI 120071: IZhang & MacFadvenll2009l) In this 
approach the range of observed times, fobs, is divided into a 
finite number (A^,) of time bins of width Afobs,/ centered on 
fobs,; (for / = l,...,Nt). That is, the ;* bin corresponds to 
fobs,/ - Afobs,//2 < fobs < fobs,/ + Afobs,//2, and there are no 
overlaps or gaps, so that fobs,/ + Afobs,//2 = fobs,/+i - Afobs,/+i/2 
for 1 < i < N, - I. For many physical systems (such as the 
ones we simulate) it is convenient to choose logarithmically 
spaced bins, with a constant Afobs,//fobs,/- If the time bins are 
sufficiently densely spaced, such that the second time deriva- 
tive (with respect to fobs) of Fy is correspondingly small, then 
Fyitoha.i, n) can be approximated by its average value within 
the /* time bin. 



(27) 



2 /-A/„b.,,,+A/„bs,,72 

^'v(fobs,/, n) = — I dtoh^Fy{tobs,n) . (30) 

Afobs,/ jA/„b,,-Af„b,.,/2 

Now given that 5[f{x - xo)] - 5{x - JCo)/l/'(JCo)l when f{x) 
has a single root at xq, we obtain 



Fy{tob,.i, n) = \ d X \ rffobs 5 f- + 

d£(z)Afobs,/ J JAf„b„-A/„b„/2 \ 1 + z c j 
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4(z)Afobs 



J J 12(1 +z 



^^^^ fd^.H 
df{z)Atobs,z.i J 



z) 

Atobs,z,i 



1 +Z 



■ + 



n ■ r 

c 

h ■ r 



Jr 



r2(l -n-j8)2 
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(31) 



r 



where H(x) is the Heaviside step function and fobs - = fobs /( 1 + 
z). 

The results of a simulation that models the dynamics of 
a physical system are naturally given at a finite number («,) 
of time steps (f.^, where j = i.e. "snapshots" of 

the dynamics. At each snapshot the values of the hydrody- 
namic variables are provided at a finite number of points, 
each at the center of a computational cell, which represents 
a finite three dimensional volume AV*^' (generally different 
from that of other cells for an AMR code). For this rea- 
son, we assign to each snapshot time t,j a finite time interval: 
(3fj.i-f.,2)/2 < t. < (fjj +f,,2)/2 and Af-j - t^^2-tz,i for ; = 1, 
(t,J_i+lj)/2 < t, < (4; + fzj+i)/2 andALj =\t,jii-t,j^i)/2 
for 2 < ;'<«(-!, and (tz.n,-i + f<-,«,)/2 < t. < (3f;,„, - fc,„,-i)/2 
and Afj,.,,, - f.„, - f,.,„,-i for j = «/. Sufficiently dense and 
well distributed snapshot times are key to the flux calcula- 
tions. Thus, the simulation provides a finite number of 4- 
dimensional space-time cells, which together cover the finite 
simulated 4-volume (the time and three dimensional volume 
covered by the simulation^). The 4-volume of the A;* 3D cell 
of the snapshot time is AVj^' - Af^jAVj^'. Given the phys- 
ical conditions in each such 4D space-time cell we can then 



calculate its local (comoving) emissivity, j'^, (under appropri- 
ate assumptions) and use equation (ISTT i in order to calculate its 
contribution to the observed flux density, Fy. The proper way 
of doing this is to calculate the fraction fij^ of its 4-volume 

Ay'.f ' that falls within each observer time bin centered on fobs /, 

1 

resulting in the following discretized version of equation (bib . 



Z7 ^, '\ (1+^)^ A 1/(4) ■^"''j'' 

(32) 

where the subscript "jk" indicates that the relevant quanti- 
ties are evaluated at the appropriate cell, centered on (f;, r) - 
{t~ j, rjit). Since the order of the summation is not important, it 
is much more convenient to evaluate the contributions of each 
4D cell according to the order at which it is stored. Since it is 
not always convenient and may cost additional computational 
time to calculate all of the coefficients fijk, one might further 
simplify equation ( l32t by attributing all of the contribution 
from any given 4D cell to a single observer time interval, cor- 
responding to that of the cell's center: 



AFyjjkih) 



(l+zf 



(4) 
jk 



v' ,jk 



dliz) Afobs,,r2^(l-,l./J^,)2 



(l+z) 



AV 



(4) 
jk 



v'Jk 



dliz) Afobs,c,, r2(l-« -4,^)2 



for 



for 



0* 



1 +z 



Afn 



n ■ rjk 



2(1 +z) 

Atobs,z,i 



(33) 



Finally, one could simplify things even further by assuming an 
isotropic emission in the fluid (comoving) rest frame, and then 
j'y,{n') - dE' /dV'dQ.'dv'dt' can be replaced by P'yJA-n where 
P'^, - dE' /dV'dv'dt' . We currently make this simplifying 
assumption. 

For 2D jet simulations, which assume an axisymmetric 
flow, the jet symmetry axis is the z-axis and it is convenient 
to choose the x-axis along the h-z plane, so that h may be 
easily expressed in terms of the viewing angle Sobs (where 
cos dobs - n ■ z). 



n = X sin Sobs + Z cos Sobs 



(34) 



Thus, in spherical (r,6,<p) or cylindrical (z,p,<p) coordinates 
(with yS^ = 0), we have 

h ■ r = r(sin 0COS (f> sin dobs + cos 6 cos dobs) = 

= p cos (f) sin ^obs + z cos 0obs , (35) 

h • ^={PrSVLlO + Pg COS 6) COS ^ sin Bobs + 

+ (/?,. COS 61 sin 61) COS 0obs (36) 
-Pp cos ^ sin 6iobs + Pz cos 6iobs ■ 



3.2. Calculation of the observed image 

The observed image can be calculated by dividing the plane 
of the sky (i.e. the plane normal to h) into bins or 2D "pix- 
els" and assigning the contribution AFyjk from each compu- 
tational 4D cell to the appropriate pixels (or pixel), where 
the conversion from flux to specific intensity (which is rel- 
evant for the image calculation) is done by using the relation 
dFy - lydS ±/d^, so that the intensity contribution to the 
pixel whose area is AS ±j would be 

^AFyjjkih) _ 



Alyjijkin) ^ dj^ 



AS^j 



(l+z)-W^' 



Jv 



v'Jk 



AS^jAtobsj T\(\-h-Pjkf 



(i+z)-^Ay]:' 



ly'jk 



(37) 



A5i./Afobs.,,; Y^.^{\-h-l3jkf 

The assignment of the contribution to the appropriate pixel 
requires a parameterization of the plane of the sky. For this 
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purpose we use a rotated reference frame denoted by a twid- 
dle, where y - y and the z-axis points to the observer (in the 
direction of n). 



X-X cos tiobs - z Sin tyobs = 

= r(sin 0COS 4> cos 0obs - cos 6 sin 6?obs) 
=p cos cos Sobs - z sin O^hs , 



y = y = r sin sin = p sin ( 



:2 + j;2 



- y sin sin 
tan0=- = 

X sin0cos(^cos0obs - cos0sin0obs 



(38) 



(39) 



(40) 



p cos COS t/obs - Z Sin Uobs, 

For an axisymmetric flow the image is invariant to y — > -y or 
equivalently to — » -0, i.e. /v(fobs, ^^y) = ^v(fobs, -^^ -y) 
and Iy{tohs,h,p,4>) - Ivitobs,n,p,-^)-^ A 2D simulation 
(whether in spherical or cylindrical coordinates) provides 2D 
snapshots of the dynamics, and each 2D computational cell 
(not counting the time dimension) needs to be transformed 
into one or more 3D cells. For the special case of an ob- 
server along the jet (or flow) symmetry axis, corresponding to 
Sobs = 0, the contribution to the observed emission (i.e. to ly 
or Fy) becomes independent of 0, which in this case is equal 
to (p, so that the image has circular symmetry {ly becomes in- 
dependent of 0) and a single bin in becomes sufficient for 
the calculation. For 0obs > 0, however, one needs to artifi- 
cially produce a large number of bins in 0, each corresponding 
to a 3D cell, which together represent a single 2D computa- 
tional region. The choice of binning should be done wisely, 
such that the Doppler factor does not vary by a large factor 
between neighboring bins (in order to calculate the observed 
radiation accurately enough) and the bin size should not be 
too coarse (as to cause excessive graininess in the calculated 
images or lightcurves), while having a reasonable number of 
bins (in order for the computational time not to be too large, 
especially for high-resolution simulations). Please note that 
since the contribution to the flux is invariant to ^ -cp, it is 
enough to choose values in the range < (p < n and give each 
resulting 3D or 4D cell a double weight when calculating Fy 
(since <pi < (p < (p2 also represents -02 < < -0i). 

3.3. Synchrotron radiation 

The main purpose of the current radiation calculations is to 
check the effect of the dynamics on the afterglow lightcurves. 
Because of this, we intentionally use a very simple model 
for the radiation mechanism (following iGranot. Piran & Saril 
Il999ah . which features synchrotron emission and ignores in- 
verse Compton scattering or its effects on the synchrotron 
emission through the additional electron cooling that it 
causes. It also ignores self-absorption, and the local emis- 
sion spectrum is approximated by a broken power-law. The 
magnetic field is assumed to hold everywhere a fraction eg of 
the proper internal energy density, e', i.e. B'^/Stt = e^e'. Just 
behind the shock all electrons are assumed to be accelerated 
into a power-law energy distribution. 



N{y<,) ^ Je'' for Je > Jm = 



p-2 

p — I I n'/UeC^ 



(41) 



^ This can also be seen from equation j36t . where the dependence on (p 
is on ly th rough cos (j>, which is invariant lo (j) ^ -(p that according to equa- 
tion (39) coiTesponds to y — > -y or (j> -(j>. 



The local emissivity P'^, is taken to be a broken power-law, 
(v'/v' v' < v' < v' 



PI 



PL 



(y'/v;)-'/^ 



y < V < y 
y' < y' < y' 



y^ < V < y, 



m 



iy' IVj^-P^I^{v' IVX^I- V > max(y;„ v^,) , 

(42) 

with the following flux normalization and break frequencies, 

512V2^ / j7 - 1 \ ql 
27 \3p-ljtn,c^' 

\2 



P' —0' 

v'.max 



{eBe')''\'^, (43) 



, 3V2;^/p-2\ 
, 27 V2^ 



efe^{e'f\n:r, (44) 



-2 



128 



0"t 



ieBe') 



/x-3/2 



(45) 



Electron cooling is treated in an approximate manner, by as- 
suming that everywhere the electrons have cooled at their cur- 
rent local cooling rate over the dynamical time, which is in 
turn approximated as f^^^ ^ r./F, so that the expression in 
equation (|45]) is simply derived from 



7c = 



3m„c 



Ao-jeBe't'^^^^ A-o-Tese't, 



\6meC 



(46) 



A more proper treatment of the electron cooling would re- 
quire following each fluid element from the point where it 
crosses the shock and the electrons are accelerated, and solv- 
ing the equation for the subsequent evolution of their energy 
distribution, accounting for their radiative losses and adiabatic 
gains or losses. This has been done analytic ally for the BMK 
self-similar solution (iGranot & Sari' 2 002[) and num ericallv 
using a ID Lagrangian code dNakar & Granotll2007h. It has 
also been im plemented in an Eulerian code (Ivan Eerten et alJ 
in a somewhat approximate fashion due to the 
difficulty in accurately tracking the electron energy distri- 
bution in each fluid element. The differences between our 
treatment of the elect ron cooling and the results presented by 
(IGranot & Sarill2()02h are shown in detail in the Appendix. 

It is also possible to use an even simpler emission model 
that ignores electron cooling altogether by assuming v', v,'„ < 

in the broken power-laws of equation (l42l) . In this paper 
electron cooling is always implemented in our calculations. 
In an accompanying paper ( iDe CoUe et alJl20Tl1) . however, in 
some cases we also use an even simpler emission model that 
ignores electron cooling altogether, 

p', { iy'ly'J" v'<v;„, 

^ = (47) 



4. APPLICATION: EVOLUTION OF A RELATIVISTIC IMPULSIVE 
BLAST WAVE 

In this section, we use our AMR+radiation code to study 
the evolution of impulsive relativistic blast waves both in one- 
dimension (ID) - a spherical blast wave propagating into 
either a uniform or a stratified medium, bridging from the 
Blandford-McKee to the Sedov-Taylor (ST) self-similar so- 
lutions - and in two dimensions (2D) - an axi-symmetric jet 
propagating into a uniform medium. 
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4.1. Self similar solution 

iBlandford & McKe^ (1 19761) studied the self-similar propa- 
gation of an ultra-relativistic spherical impulsive blast wave 
in a medium with a density 



Pk(r) 



(48) 



They showed that an appropriate choice of the similarity 
variable is 



X ^1+2(4- k)r- 



(49) 



where r and R are the radial position (in polar coordinates) 
of the fluid element and of the shock front respectively; Fsh is 
the Lorentz factor of the shock front, which as that of the fluid 
(and all of the velocities) here is measured in the rest frame 
of the upstream medium ahead of the shock, and it is related 
to the the Lorentz factor of the shocked fluid just behind the 

shock front by TOr = 1) = Tsh/ V2. IBlandford & McKed 
(1 19761) showed that the position of the shock front is given by 



R^ct 



1 



1 



2(4 - k)r^ 



and its Lorentz factor can be written as 

(17 -4-k)E 
^7Tpt{R)c^t^ 



^ sh - 



(50) 



(51) 



where pk{R) = A^R^^ is the density of the ambient (un- 
shocked) medium at the position of the shock front, and E 
is the energy in the blast wave. 

The lab frame time corresponding to a given Lorentz factor 
of the shock front is therefore (see equations ( |48] |. (ISOl l and 
dSTl i) given by 



^ ^ 1 
c c 



(11 -4k)E 



llO-k) 



(52) 



The post-shock Lorentz factor F, proper rest-mass density 
p, and pressure p, are given by 



V2 

p^2'I^Pk{R)V,^ 



-(10-3<:)/[2(4-i)] 



(17-4*)/[3(4-i:)] 



(53) 

(54) 
(55) 



The relativistic blast wave typically begins to slow down 
when it sweeps up an amount of mass with a rest-mass en- 
ergy of order of the kinetic energy of the blast wave. That 
corresponds to a distance (Sedov length) of 

l/(3-A-) 



(3 - k)E 



(56) 



where the jet energy E is the energy (excluding rest energy) 
in the flow. For a non-spherical flow, or a jet, to zeroth 
order E in equations ( BTT i and ( |52] | can be replaced by the 
local value of the isotropic equivalent energy in the flow, 
fiiso = Ajx{dE IdQ), as long as it does not vary significantly 
over an angular scale of the order of the inverse of the lo- 
cal value of the Lorentz factor of the fluid just behind the 
shock. In particular, for a double-sided conical wedge of half- 
opening angle 9q taken out of the BMK solution (or a uni- 
form sharp-edged jet), which we later use as the initial condi- 
tions of our 2D simulations, we have E - {I - cos 0o)£iso ~ 



(0Q/2)£'iso ~ 2 X 10^'(6'o/0.2)^£'iso,53 erg, where we have used 
a fiducial value of Ei^o - lO^^iSiso.ss erg, typical for long du- 
ration GRBs. Whether it is more appropriate to use E or E^^o 
in equation ( |56] | for such a jet, i.e. at which distance from the 
origin it becomes Newtonian, is a non-t rivial question, whic h 
is addressed in an accompanying paper (iDe Colleet al.l2011h . 

In the non-relativistic limit, the self-similar b ehavior of the 
blast wave is described by the Sedov-Taylor dSedovl [T959I: 
iTivlor 1950) self-similar solution, with the position of the 
shock wave given by 



R 



Ak 



(57) 



and the shock velocity given by v,/, = dR/dt cc f-(3-*)/(5-*) 
Approximated expressions for the post-shock density, pres- 
sure and velocity profiles in the ST regime are given, e.g., by 
[Petruk (2000). As there are not analytical solutions for the 
scaling of density, pressure and velocity in the post-shock re- 
gion, it is not possible to find a simple analytical expression 
for a^. Based on the simulations presented in ^4.2t , we find 



a 



l/(5-k) 



= 1.15, 1.04, 0.78 for /t = 0, 1 , 2 respectively. 
4.2. Initial conditions 



In this paper, we perform a series of ID (with ^ = 0, 1,2) 
and 2D (with ^ = 0) simulations of the propagation of impul- 
sive blast waves, including the transition from the relativistic 
to the non-relativistic phase. All simulations employ spheri- 
cal (polar) coordinates, and using the HLL method (see 32.2| i 
for the flux calculation. The multi-dimensional simulations 
for the cases k - 1, 2 are presented in an upcoming paper 
jDe CoUe et aLllMl . 

The initial conditions of the problem depend on the values 
of the following parameters: the isotropic energy of the blast 
wave, Ei^o, the initial Lorentz factor of the jet shock front, 
Fsh,o, the density profile of the external medium (that is, the 
values of k and of the normalization factor A^) and the jet 
initial half-opening angle, 0o (in the 2D case). In all the sim- 
ulations, the initial profiles of density, pressure and Lorentz 
factor (radial velocity) in the post-shock region are set from 
the BMK self-similar solutions, given by equations (I53b-(l55b. 
We initialize the density of the ambient medium (in the case 

= 0) as Ao = po = no'«p = 1.67 x 10"^"* g cm"^*, and the 
pressure as p - tjpqc^, with rj - lO"'". The value of rj does 
not affect the outcome of the simulation as long as the Mach 
number remains large, i.e. Ai ~ 77"'^~v'sh/c » 1. As the sim- 
ulation continues to evolve well into the Newtonian regime, 
this condition corresponds to Vsh » 3 (?7/10"'°)'^^ km s"'. 

In a first set of simulations, we study the deceleration of 
mildly relativistic impulsive blast waves bridging from the 
BMK to the ST self-similar solution. In th e case ^ = 0, the 
initial conditions are similar to those used bv lvan Eerten et all 
( I2OIOI) . To determine the density profile in the cases k - 1,2, 
we fix the Sedov length (equationl56b as Ls{k) = Ls(k = 0): 



is 



(3-k)E 


l/O-k) 


3E 


1/3 


4nAkC^ 




AnAoc"^ 





(58) 



and derive an expression for A^ as A^ = AoLf (3 -A:)/3. There- 
fore 

3-k( rV'' 
P = A„^(-) . (59) 

We further assume Ei^o = 10^^ ergs, corresponding to a Se- 
dov length ofLj = 1.17 x 10'^ cm, and a Lorentz factor of 
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the shock of Fsh.o = 10. To properly cover the deceleration to 
non-relativistic speeds (especially for the case k = 2), we use 
a large spherical box of radial size = 3 x 10^'^ cm (corre- 
sponding to a size of ^ 256Ls)- The simulation is stopped at 
ffin - 500 yrs. 

In the case k - 0, the simulations begins at fo = 1.19 x 
10^ s, with a jet shock located at Rq - 3.56 x 10'^ cm. 
The case k - I corresponds to an initial time and jet 
shock radius given by fi/fo = Ri/Ru = 0.53. The case 
k = 2, corresponding to a steady spherically symmetric 
wind, has 12/ to = R2/R0 - 0.074. The values assumed 
for the spherical wind can be compared with those ob- 
served for Wolf-Rayet stars, which winds have large mass- 
loss rates of M a; lO"'' - 1 0'^ Mpyr ' and velociti es v'„ x 
1000 - 2500 km s"' (e.g., IChiosi & Maededll986h . giving 
n^(r) « 0.45(r/10"^cmr2(M„,./3 X 10"^ Moyr-i)(v„r/2 x 
lO-'kms"') cm""*, which is very similar to the one used in 
the simulations. 

The AMR code uses a basic grid of 1000 cells with a maxi- 
mum of 18 levels of refinement, corresponding to a maximum 
resolution of Arn,in = 2.3 x 10'^ cm. In a uniform grid code, 
the same resolution would be achieved by using 1.3x10*^ cells. 

In a second set of simulations, we test the radiation code 
by running simulations of highly relativistic decelerating blast 
waves (limited to the case k = 0) both in ID and 2D. In 
these simulations, we assume an isotropic energy of ^iso = 
10^^ ergs, corresponding to a Sedov length of = 2.51 x 
10 cm, and a Lorentz factor of Tsh.o = V2 X 20. The sim- 
ulations begins at to - 1.277 x 10^ s, with the shock initially 
located at Ro - 3.83 x 10'^ cm, and ends at f^n - 150 yrs. To 
properly study its lateral expansion, an initial opening angle 
of 60 = 0.2 rad (in the 2D case) is assumed for the jet. 

The spherical box has a radial size of Lr - 1.1 x lO'** cm 
and angular size (in the 2D simulation) Lg = n/l. The AMR 
code uses a basic grid of 100 cells along the radial direction 
both in ID and 2D, and 4 cells along the 6 direction in the 2D 
simulations. We run a series of simulations varying the maxi- 
mum number of refinement levels. The lowest resolution sim- 
ulation uses 10 maximum levels, while the highest employs 
18 levels of refinement in ID and 15 in 2D, corresponding to 
a maximum resolution of Armin = 2.1 x 10" cm in ID and 
Ar^in = 6.7 X 10'^ cm, A^min = 2.4 x 10"^ rad (along r and 6) 
in 2D. The structure of the grid at the beginning of the simu- 
lation is shown in Figure |2] for the 2D run. In a uniform grid 
code, the same resolution would be achieved by employing 
5.2 X 10^ cells in ID, and ~ 10" cells in 2D. 

To keep approximately constant the resolution of the rela- 
tivistic thin shell A oc the maximum number of lev els 
of refinement Meveis is decreasing with time (lGranoi2007h as 
Meveis = max[7,Mevds.o - (4 - fc)log(f/fo)/log(2)]. We re- 
fine our adaptive mesh based on rest mass density and energy 
gradients. The ID simulations run in at most a few hours on 
a normal workstation, while the 2D simulations need a few 
days on ~ 100 processors. 

4.3. One-dimensional simulations of trans-relativistic blast 
waves propagating in a stratified medium (k — 0, 1, 2j 

During its deceleration, the shock front is typically re- 
solved with 3-4 cells (Figure[3]l, as is the case for most mod- 
ern Eulerian shock capturing schemes. The normalized lab 
frame density behind the shock, given from the relativistic 
Rankine-Hugoniot conditions for strong shocks: plpk{R)T - 
(y + l/r)/(y - 1), remains approximately constant during 




0-02 o.oi 0.06 O.oa 0.10 

R (xlO'ia cm) 



Fig. 2. — Adaptive grid structure for tlie initial condition of the two- 
dimensional simulation of a relativistic blast wave. The blue color indicates 
the post-shock region, while the green area represent the ambient (unshocked) 
medium. 




r [cm] 



Fig. 3. — Density profiles normalized to the Lorentz factor and the local 
value of the ambient medium density. The curves shown in the upper panel 
(k = 0) correspond to r = r,, = 137 days, and r = 1,3, 20, 100, 500 yrs. The 
central panel {k = i) includes also the profile at ?o = 81 days. The bottom 
panel (k = 2) includes also the profile at fo = 20.3 days. The horizontal red 
line indicates plpk(r)T = 4. 

the t ransition from relativistic to non-relativistic regime s (see, 
e.g.. lBeloborodov & Uhmll2006l: Ivan Eerten et al.ll2010l) . Fig- 
ure [3] shows that in fact plpi.(R)T ^ A ?A different times and 
for different values of k. The drop of the density profile in 
the post-shock region approximately follows the BMK self- 
similar solution, and is therefore less steep with larger k (see 
equationl54li. This figure also shows that the deceleration pro- 
cess is slower for a more stratified medium. 

Figure |4] shows the evolution of the shock front radius 
for different density stratifications. Both the ultra-relativistic 
(with /?sh ~ ct) and the non-relativistic (R oc ^2/(5-*:)-^ ^jj^jy]-. 
ical self-similar solutions are properly recovered by the sim- 
ulations. As shown e.g. by Ivan Eerten et al.l (1201 Ol) for the 
case A: = 0, the transition from relativistic to non-relativistic 
phase happens on scales much larger than L,. If for instance 
we estimate from Figure |4] the time it takes for the relativistic 
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Fig. 4. — Position of the sliock front for tlie three cases/: = 0, 1,2 
(up to bottom) along with the ultra-relativistic (/?,;, = ct) and the Sedov- 

Taylor (Rji, = ^fftiJisoJ^/At.^ ) regimes. The Sedov-Taylor curves as- 
sume = 1.15, 1.04, 0.78 for i = 0, 1, 2 respectively. The gray curves 
are computed from a semi-analytical approximation based on energy conser- 
vation (see the text for a detailed description). 
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Fig. 5. — Number density profile (in the lab frame) for dift'erent resolutions 
for the case k = at the beginning of the simulation (f = 148 days , upper 
panel) and ? = 156 days (bottom panel). 



blast wave to slow down to non-relativistic speeds based on 
the intersection between the relativistic and non-relativistic 
self-similar curves, we obtain values of ~ 0.9 x 10^ days, 
1.2 X lO-' days and 2.7 x lO-' days for A: = 0, 1, 2 respectively. 
These values are much larger than those computed by using 
the Sedov length (Piran 2005) Jnr ~ Lnr/c - 450 days. 

This result, together with the scaling of position, Lorentz 
factor and the shock velocity as a function of time and shock 
radius, can be easily understood by a simple analytical argu- 
ment involving the conservation of energy. In fact, the energy 
is given in the ultra-relativistic regime by 

^ = T^^^^'^'^'ry (60) 
17 - 4k 

and in the non-relativistic limit by 

E = ^IzJl^Akv'-R'-' (61) 

As the energy has a common scaling in relation with the 
other physical parameters iv,R), differing only in the constant 
of proportionality, a simple interpolation between the two lim- 
its is given by 

, = „»^r%.=(^,'.(5^(i-,?,) (62, 

This equation can be easily written as function of velocity 
as: 



1 + CMRiR/L^y-' + V[i - c^R{R/uy-']- + 4cR{R/L,y-k ' 

(63) 

where cr = and cnr = ^^-^i^- This expression ap- 

proximately gives the dependence of v (or F) on the shock 
position, for every choice of the blast wave energy and den- 
sity stratification. For instance, at R ~ Ls, equation ( l63T l 
gives Vsh ~ 0.83,0.85,0.89c (or M = FyS ~ 1.46,1.64,1.99) 
for k - 0,1,2 respectively. At this radius (and time) the 
shock is therefore still relativistic, and the ST solution is not 
valid. The exact determination of depends, however, on 
the definition of the trans ition between the relativistic and th e 
non-relativistic flow (e.g. lRamirez-Ruiz & MacFadvenl20Tol) . 



If for instance we define Jnr as the time where the asymp- 
totic BMK solution and the ST power laws are equal (i.e. 
ct/L, = [4TOj(3-/t)] !/"-*>), wegetf ~ 2.1fNR ~ 9x10' days 
(k = 0), f ~ 3fNR ~ 3.4 yrs {k = 1) and t ~ 6fNR ~ 7.5 yrs 
(k = 2). At this time the blast wave is nonetheless still mildly 
relativistic (fi - 0.51,0.56,0.63) and the ST solution is not 
valid. If on the other hand we assume that the ST solu- 
tion becomes valid at a fixed (somehow arbitrary) speed of 
v/c < 1/3, we get t ~ 3.6rNR ~ 1-6 x 10^ days (k = 0), 
t ~ 7fNR ~ 8.6 yrs (k = 1) and t ~ 48fNR ~ 59 yrs (k = 2) 
(Figure |4]i. 

Equation (|63l l. when rewritten in the form dR/dt - /3{R), 
admits a complex solution t = f{R) in terms of Appell hy- 
pergeometric functions. The time dependence of the shock 
position R - R{t) has been therefore more easily derived by 
numerical integrating equation ( |63] l. and it approximates the 
position of the shock computed from the numerical simula- 
tion within a maximum difl'erence of 1, 2, 5% (for A: = 0, 1, 2 
respectively). 

While Figures [3] and |4] clearly show the validity of our 
implementation for mildly-relativistic and non-relativistic 
speeds, reproducing the correct BMK self-similar scaling dur- 
ing the early stages of the simulation, when F > 10, is much 
more challenging. 

Figure|5]shows the initial density profile (for the case ^ = 0) 
in the region around the position of the shock. A very large 
number of levels of refinement must be used to properly ini- 
tialize the density, pressure and Lorentz factor in the post- 
shock region. For instance (Figure |5] upper panel), the initial 
steep density profile is recovered with errors less than 10% 
only by using resolutions corresponding to > 18 levels of re- 
finement. 

While the BMK self-similar solution represents an exact so- 
lution of the SRHD equations in the ultra-relativistic limit, 
the particular discretization employed may not be the ex- 
act (numerical) solution of the discretized equations. As a 
consequence, the relaxation towards the numerical solution 
passes through the development (see Figure |5] bottom panel) 
of a spurious numerical "precursor" propagating in front of 
the BMK shock if insufficient resolution is used. While the 
size of the precursor shock drops effectively with resolution, 
it also produces a quick drop in the maximum Lorentz fac- 
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Fig. 6. — Maximum Lorentz factor in the post-shock region (measured in 
the lab frame) as a function of time. The simulations (with k = 0) start at 
t ~ 147 days with a Lorentz factor of 20. The curves shown correspond to 
the expected BMK self-similar solution, 14, 16, 18 lev els of re finement with 
a fixed maximum level of refinement, and results of the lZhang^ife MacFadyen 
120091) two-dimensional simulations. For 14 levels of refinement we show 
two curves, either with (a) or without (b) decreasing the maximum level of 
refinement with time. Each label on the x-axis corresponds to the time when 
the maximum resolution drops by a factor of two, so that for instance the sim- 
ulation with initially 14 levels of refinement drops to 13 levels after 175 days, 
12 after 209 days, and so on. 



tor behind the shock (due to the spreading of the initial F 
peak, see Figure |6]l. The Lorentz factor eventually converges 
to the correct BMK solution at F ~ 10 at the largest reso- 
lution used (18 levels of refinement). Figure |6] also shows 
the effect of decreasing the maximum level of resolution 
during the evolution of the simulation dGranot et al.l 120011: 
IZhang & MacFadyeri '2009': D e Colle et al.ll201 Ih . As can be 
appreciated from Figure |6] the decrease in the resolution pro- 
duces a slower conv ergence to the BMK solution . The time 
evolution of F from IZhang & MacFadvenI ( |2009|) . included 
in Figure |6] is similar to our low resolution (14 levels) one- 
dimensional simulation, corresponding approximately to the 
resolution ac hievable in multi-dim e nsiona l simulations. The 
noise m the IZhang & MacFadvenI (l2009l) curve is due to a 
larger temporal sampling. A proper treatment of the tiny ultra- 
relativistic post-shock region would require a larger resolution 
or alternatively a much less diffusive method as e.g. high or- 
der (cou pled to high resolution) or Lagra ngian-Eulerian meth- 
ods (e.g. lKobavashi. Piran & Sarill 19991) . 

The specific numerical resolution required is determined by 
the relevant structure one needs to resolve. The hardest to 
resolve, in our case, is the initial BMK shell (A) at the initial 
time (fo) or radius (Ro)- Its effective width does not have a 
unique definition, but it can be parameterized as 



Aa - a- 



Ro 



r?h(^o) 



(64) 



where the numerical factor a can be evaluated using the BMK 
self similar solution Defining Ao as the width of the region 
behind the shock that contains a fraction / of the total energy 
(E) or rest mass (M), respectively, results in 



2(4 -k) 



as = 



3(4-^) 
17-4;t 



4-k 
3-k 



(65) 



' Note that if one uses r{Ro,x = 1) instead of rsii(Ro) in equation (64) 
then the value of the numerical coefficient a would be smaller by a factor of 





t [days] 

Fig. 7. — Light curves (at lo'''^''^ Hz from up to bottom) from simula- 
tions at different resolution, either including (right panels) or not including 
(left panels) the contribution from the synthetic lightcurve emitted from a 
Blandford-McKee self-similar blast wave with Lorentz factor between 200 
and 20. The synthetic light curve (labeled BMK in the Figure) emitted from 
a Lorentz factor between 1 and 20 (or between 1 and 200) is also shown in 
the left (right) panels of the figure. 



For / = 1 /2, this gives fl£ = 0.0789, 0.103, 0.147 and am = 
0.190, 0.305, 0.750 for = 0, 1,2. 

One can then similarly express the numerical resolution in 
terms of a parameter 0^5, 



^0 



ri(^o) 



(66) 



where Armin is the smallest resolution element in the radial 
direction. Previous 2D jet numerical si mulations wi t h sim- 
ilar initial conditions used ^ = 0. In iGranot et al.l (12001 1) 
the initial resolution was ra ther poor, Ares = 0.69, while in 
IZhang & MacFadvenI (120091) it was significantly improved, 
cjres = 0.12. Here we use Ares - 0.014 for k - Q, which 
represents an order of magnitude improvement. For k = I and 
2 we have fl,e,s = 0.022 and 0.087, respectively. 

Figure [T] shows that the light curve, computed by post- 
processing the results of the simulations with our radiation 
code, converges quickly except for fobs S 0.5 days, where 
part of the flux, which should be generated from regions with 
r ~ 20, is shifted to a lower fobs- That can be in part com- 
pensated by adding the contribution coming from the jet de- 
celerating with 20 < r < 200, computed by mapping in the 
radiation code a BMK self-similar solution. As shown in Fig- 
ure |7] (right), the sum of the synthetic flux with 20 < F < 200 
and the flux computed from the results of the simulation with 
1 < F < 20 produces a valley (shallower for increasing res- 
olutions) for fobs ~ 1 day. This artificial feature is due to 
relaxation from the initial conditions to the numerical solu- 
tion, and gradually disappears as the resolution is increased. 
A comparison between the lightcurve computed from the ID 
simu lation (with k = 0) an d the semi-analytical calculations 
from lGranot & Saril (120021) is shown in the Appendix. 

4.4. Two dimensional simulations for k — 

Figure [8] shows snapshots representing the early evolution- 
ary stages of the jet density. During the relativistic phase, 
there is only modest lateral expansion. As portions of the 
jet expand laterally, a rarefaction front moves towards the 
jet axis. The strong shear present at the contact discontinu- 
ity drives shearing instabilities that have however a negligible 
effect on the shock dynamics and afterglow radiation com- 
ing from the jet. At the jet break time f = fjB ~ 8.7 yr. 
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Fig. 8. — Lab-frame density stratification snapsliots of the 2D simulation at 147 days (left), 256 days (center), 372 days (right panel). 



the lateral expansion becomes more vigorous, and at later 
stages (on times » Jnr) the jet slowly converges to a spher- 
ical shape. Although it is not possible to make a quantita- 
tive comparison, our results qualitatively resemble those of 
IZhang & MacFadvenI (120091 see their Figure 2 for a direct 

comparison), as well as those of iGranot et alj (12001b. 

Wh ile theoretical arguments dGruzinovl 120001 : IWang et al.l 
I2OO2I) seems to indicate that the shock front should be stable 
to linear perturbations for either a uniform or a wind den- 
sity profile of the ambien t medium, recent simulations by 
iMeliani & KeppensI (1201 Oh observe the development of insta- 
bilities in the shock front. The deve lopment of similar insta- 
bilities is also observed by iDe Coll e et al. (2011) relative to 
the case of a stratified medium with k - 2, while it is not 
observed in the simulations presented i n this paper (despite 
using the same HLL Riemann solver as iMehani & Kepperisl 
I2OIOI and similar initial conditions), consistently with the re- 
sults by Zhang & MacFadyen (2009). The different results in 
the simulation seems t o imply a numerical origin for the in- 
stabilities observed by iMeliani & KeppensI (1201 Ol) . although 
further investigation is needed to better understand the prob- 
lem. 

The afterglow light curves computed from our 2D jet sim- 
ulation assume that the observer is located along the jet sym- 
m etry axis (0obs - 0). To facilitate comparison with the results 
of IZhang & MacFadvenI (l2009h . we choose the same param- 
eters for the afterglow calculation: eg = = 0.1, z = 1 and 
p - 2.5, in addition to the same values for the parameters to 
determine the hydrodynamics (Zsiso = 10^^ erg, next = 1 cm"^* 
and 0() -0.2 rad). 

As in the ID case, the afterglow emission (Figure |9]l shows 
a shallow valley at f < 1 day, due to a lack of resolu- 
tion into the region immediately behind the high relativistic 
shock. Figure |9] (bottom panel) shows a comparison with a 
2D "wedge" (computed by using a ID simulation mapped on 
a wedge with 6 < 0.2; the finite resolution of this ID simula- 
tion is affecting the lightcurves at the earliest times as shown 
in Figure |9]l. Before the jet break time, the 2D light curve 
from the simulation is very similar to that from a 2D wedge 




Fig. 9. — Afterglow emission at lO''"''^*'" Hz from the 2D simu- 
lation compared with a 2D wedge (bottom panel) and the results from 
IZhang & MacFadvenI j2009l) (upper panel). 



with the same (initial) isotropic energy, indicating that little 
sideways expansion takes place befor e the jet break, in agree- 
ment wi th previous analyti cal (e.g.. lRhoadslll999l) and nu- 
merical (iGranot et al.ll200l"l) results. After the jet break time, 
however, the flux from the 2D simulation becomes lower than 
that for the corresponding wedge, and the difference between 
the two gradually increases with time, as the lateral spread- 
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Fig. 10.— Spectra at = 0.1, 1, 10, 100, 1000 days (black, red, green, 
blue, puiple) 



ing of the jet gradually increases during the relativistic phase 
and then more rapidly during the Newtonian phase (until at 
very late times spherical symmetry is approached). Our cal- 
culated afterglow emission and spectra agree very well with 
IZhang & MacFadvenI (120091) (Figure |9] upper panel and Fig- 
ure [Tol l both in the flux before and after the jet break. 

5. CONCLUSIONS 

In this paper, we have presented a detailed description of 
the new state-of-the-art adaptive mesh refinement, relativistic 
hydrodynamics code Mezcal-SRHD and of the radiation code 
used to compute the synchrotron emission from the output of 
the hydrodynamics simulation. The proper implementation 
of the SRHD algorithm has been verified by running stan- 
dard one- and multi-dimensional tests which are presented in 
the Appendix. The code has been applied to the study of the 
propagation of ultra-relativistic impulsive blast waves both in 
one and two-dimensional spherical coordinates. 



We have studied for the first time the deceleration of rel- 
ativistic impulsive blast-waves in one dimension propagat- 
ing in a stratified medium and find that the deceleration to 
non-relativistic speeds happen on scales /?nr from a few (for 

= 0) to several times larger than the Sedov length L^. Tak- 
ing /?NR as the radius where Rsjit) = ct gives the expres- 
sion /?nr/As = [^^o'k/{3 - k)]^^^-^~^\ which iflustrates how 
7?NR/is increases with the degree of stratification of the am- 
bient medium where the shock is propagating. These results 
have been described in detail using a simple semi-analytical 
formula, derived from energy conservation, which gives the 
correct scaling of the position and velocity of the shock as a 
function of time. 

The results obtained by the radiation code were validated 
by a comparison with semi-analytical results, and with those 
obtained in previous numerical works. We have also shown 
that while the resolution is a key factor to properly recover the 
correct dynamical evolution of the system (with some of the 
parameters not yet converging, e.g. the shock Lorentz factor), 
when the contribution from the radiation produced by the jet 
before the onset of the simulation (in our case 20 < Fsh/ V2 < 
200) is included in the calculation, the resulting light curve 
becomes much less sensitive to the exact resolution. 

In an upcoming paper, we will extend the results of the sim- 
ulations presented here to include multi-dimensional simula- 
tions in a stratified medium. The study of the contribution of 
the magnetic field on the jet dynamics and afterglow radiation 
is left for future works. 
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APPENDIX 

EVALUATING THE APPROXIMATIONS USED IN THE ELECTRON COOLING FREQUENCY ESTIMATION 

A comparison between the lightcurve computed by ma pping in the radiation code a blast wave described b y a BMK self- 
simila r solution and the semi-analytical calculations from lGranot & S aril (l2002h is shown in Figure lAll While iGranot & S aril 
(I2OO2I) obtained smooth spectral breaks, for simplicity we use here their broken power-law prescription (without synchrotron 
self-absorption). In that work the afterglow emission from the BMK solution is calculated for an exact local synchrotron spectral 
emissivity while analytically calculating the electron energy distribution everywhere by following its evolution from the shock 
front (where it is assumed to be a pure power-law) due to radiative and and adiabatic cooling. The light curve computed by using 
a simplified emission model (equation l47Ti that neglects electron cooling altogether is an very good agreement with the GS02 
semi-analytical results (see Figure IaTI) . The light curve computed by using an approximated electron cooling presents three 
breaks at low frequencies (corresponding to the transitions v < < v,„ with the scaling Fy oc f '/^ ^ v < v,„ < Vc with Fy oc f'^^ 
— > v,„ < V < Vc with Fy oc f3(i-p)/4 _^ y^^^ < Vf < V with Fy oc f(2-3p)/4-) breaks at high frequencies (corresponding to 

V < Vc < v,„ with Fy oc f'/^ — > Vc < V < v,„ with Fy oc r^^'^ — > v^. < Vm < v with Fy oc f(2-3p)/4-)^ ^^j^ noticed in Figure 
lAll our estimation of the cooling break frequency Vc assuming that the electrons cool at their current local cooling rate over the 
dynamical time (see equations l42li underestimate the cooling frequency determined by GS02. For instance, an increase in Vc of 
a factor of 4 produces a better agreement with the GS02 results (Figure [ATI right panel). It is worthwhile to stress that, while 
the mapped BMK light curve and the GS02 results are applicable only for (highly) relativistic flows, the light curve computed 
from the numerical simulations is valid during all the deceleration of the flow to non-relativistic speeds. Finally, we notice that at 

V < 10'' Hz, self-absorption dominates and the light curves computed with our simple prescription are inaccurate. 

NUMERICAL TESTS 

We present in this section a series of one-dimensional shock tubes and multi-dimensional tests. 
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Fig. A1. — Comparison between light curves (at IQ^/' '/'-^/'^ /'^ Hz) computed from a Blandford-McKee self-similar blast wave with Lorentz factor between 
600 and 1, and the semi-analytical results from (Granot & Sari: 2QQ3) - Left panel: Simple emissi on m odel excluding electron cooling (equation |47) . Center. 
Light curve computed by using an approximated emission model for the electron cooling (equations|42). Right: The same as the center panel, but with a cooling 
frequency four times larger. 



One-dimensional shock tubes 

Shock tube tests are used as standard tests as they are simple to implement and the exact analytical solution is known. The 
tests were performed using a grid with size < x < 1, with an initial discontinuity at x = 0.5. Here and in the following, we 
refer to the left/right hand side of the discontinuity with the suffixes £/r. In all the tests, we use a grid with 50 cells at the lowest 
level, with 4 levels of refinement, corresponding to an eff'ective resolution of 400 cells. We also make high resolution runs of the 
same tests, employing 400 cells at the lowest level, with 4 levels of refinement, corresponding to an effective resolution of 3200 
cells. The Courant number if fixed equal to 0.8 in all tests, with a final integration time of f = 0.4. The politropic index is fixed 
equal to 4/3 in the first shock tube test and 5/3 in all others tests. As described in the following, in all the tests the exact solution 
is properly recovered. 

The first test consist of a low-relativistic flow with a left state given hy pi - \, pi - \, vi - 0.9, corresponding to a Lorentz 
factor of F »; 2.3, and a right state given by pr - 10, pr - \,vr — 0. The evolution of this shock tube co nsists of two shocks 
and a statio nary contact discontinuity. Small oscillations, similar to those observed by previous authors (e.g. lLucas-Serrano et alJ 
|2P04; Wan g et ani2008h . are present in the post-shock region. 

The second shock tube consists of a low-relativistic flow with a left state given by pi - 10, pz, = 1, v^, - -0.6 and a right 
state given by pR - 20, pr - 10, vr - 0.5. In this test, two rarefaction wave are produced, together with a left moving contact 
discontinuity. Both rarefaction wav es are properly recovered, while the contact discontinuity is smeared over ~ 10 cells. 

The last two tests are taken from lDoriaa ( Il998l) . and refer to blast wave explosions. The third shock tube consists of of a left 
state given with pi - 40/3, pi - 10, and a right state given by pR - 10"^, pr - 1, while in the last test the left state is given by 
Pl - 1000, pi-\, and the a right state is given by pR - 0.01, p« = 1, The large pressure gradient produces a mildly relativistic 
shock (test 3) and a highly relativistic shock (test 4) with F » 6. As can be seen in Figure IB II the solution consists in both 
cases of a strong shock moving to the right and a rarefaction wave moving to the left. No oscillations are present in the solution. 
The shock is resolved within ~ 4 cells, while the contact discontinuity is smeared over several cells. That is expected, due to 
the intrinsic diffusive properties of the HLL schemes. In the second blast wave problem, the size of the thin dense shell in the 
post-shock region consists of only x 4 cells with the resolution employed. As a consequence, the exact value of the density is not 
recovered at low resolution. However, this region is properly resolved in the high resolution run. 

Multi-dimensional tests 
Relativistic 2D Riemann problem 
This test has been st u died in the non-relativistic case by iLax & Liul (Il998h . and extended to the SRHD case by 



Del Za nna & Bucciantini' (20 02|). It has been widely used recently as a test for multi-dimensional SRHD codes (e.g 
Lucas-Serrano et al.n2004^ iWang et al.ll2008l) . The computational domain (at f = 0) is divided in four regions 



{p,v^,vy,pf"^ = (0.1,0,0,0.01) 



(P, V^, Vy, P) 

(p,V_„\'y,p) 

(p, V^, Vv, p)' 



SE 



= (0.1,0.99,0,1) 
^'^ = (0.5,0,0,1) 
(0.1,0,0.99,1) 



if 
if 
if 
if 



X > 0.5,y > 0.5 
X < 0.5,y > 0.5 
x>0.5,y<0.5 
X < 0.5,3/ < 0.5 



We use a uniform grid with 400 x 400 cells, an adiabatic equation of state with constant y - 5/3, and outflows boundary 
conditions. The simulations ends at f = 4. To better resolve the contact discontinuity, a more compressive MC limiter is used 
here. The results are shown in Figure lB3] The initial discontinuities across the four regions of the grid produce stationary contact 
discontinuity (with jumps in transverse velocities) between SW-NW and SE-SW, and shocks between NE-NW and SE-SW. These 
shocks produce an elongated jet-like structure on the diagonal. These features, together with the curved shock in the SW region, 
are qualitatively similar to those obtained by previous authors. 
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Fig. B2. — Logarithm of the density for the Relativistic 2D Riemann problem att = 0.4. Thirty equally spaced contours are plotted in the Figure. 
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Fig. B3. — Logaiithm of the density for the Emei'y step problem at / = 4.26. Thirty equally spaced contours are plotted in the Figure. 



Emery step 

The "Emery step" test has become a standard test both for non-relativistic and relativistic hydrodynamics codes, and it consists 
of a wind moving through a tunnel. Our initial conditions closely follow those by Lucas-Serrano et al. (2004). A relativistic flow 
moves initially horizontally with velocity - 0.999c, corresponding to a Lorentz factor of F s; 7. The density is initially fixed at 
p - \ A everywhere, with a pressure of p - 1/9 and an adiabatic index of y = 7/4, corresponding to a Newtonian Mach number 
of M = 3. The size of the tunnel is < x < 3 and < y < 1. A step is located in the region defined by x > 0.6, y < 0.2. 
Inflow boundary conditions (with the same values used to fill the tunnel initially) are fixed at the left boundary. Outflow boundary 
conditions are fixed at the right boundary, while reflecting boundary conditions are fixed at the upper, lower and step boundaries. 
We use a uniform grid with 240 x 80 cells, with the HLL method coupled to the MC limiter. 

Figure |B3] shows the density stratification at f = 4.26. As the relativistic flow collides with the step, a reverse shock is formed. 
This sh ock front is reflected from t he upper boundary forming a stationary Mach stem. The results of this tests are similar to 
those of lLucas-Serrano et al.l (l2004h . 
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